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Abstract. 

Dense suspensions of small strongly interacting particles are complex systems, which are rarely 
understood on the microscopic level. We investigate properties of dense suspensions and sediments 
of small spherical AI2O3 particles in a shear cell by means of a combined Molecular Dynamics 
(MD) and Stochastic Rotation Dynamics (SRD) simulation. We study structuring effects and the 
dependence of the suspension's viscosity on the shear rate and shear thinning for systems of varying 
salt concentration and pH value. To show the agreement of our results to experimental data, the 
relation between bulk pH value and surface charge of spherical colloidal particles is modeled by 
Debye-Hiickel theory in conjunction with a 2pK charge regulation model. 



PACS numbers: 82.70.-y, 47.ll.-j, 47.57.Qk, 77.84.Nh 
I. INTRODUCTION 

We simulate colloids of silt particles, for 
which in many cases the attractive Van-der- 
Waals forces are relevant. These colloids 
are sometimes called "peloids" (Greek: clay- 
like). In contrast to clays consisting of thin 
plateletspj our particles are in first approx- 
imation spherical particles. For real clays 
the particle shape and their orientation is of 
relevance 0, El For silt particles on 
the other hand the description is less complex. 
However, due to the particle size of micro me- 
ters and below, the interplay of diffusion, elec- 
trostatic repulsion, van der Waals attraction, 
and hydrodynamics still renders the suspension 
a very complex system. Colloid science tries 
to investigate the properties of such suspen- 
sions and there is a vast amount of literature 
on this subject 0, H H 03 OH Colloids in 
general have various applications ranging from 
food industry over paintings and cosmetic prod- 
ucts to applications in photographic processes. 
Particles with well defined properties can be 
used to investigate general properties of soft 
condensed matter like gelation or crystalliza- 
tion on a larger length scale than on the atomic 
level. Especially attractive interactions (deple- 
tion forces as well as van der Waals attrac- 
tion) have drawn the attention in the recent 
years[H 030003. 

In soil mechanics real 



samples, e.g., of sediments, can be less charac- 
terized and therefore it is more difficult to gain 
a microscopic picture from which general prop- 
erties can be derived. Therefore we have chosen 
a synthetic AI2O3 powder suspended in water 
as a model system for silt. The particle diame- 
ter is 0.37 //m. 

AI2O3 is not only a cheap testing material 
for investigations related to soil mechanics, but 
it is also an important material for ceramics. 
In process engineering one of the basic ques- 
tions is, how to obtain components of a prede- 
fined shape. Wet processing of suspensions, fol- 
lowed by a sinter process is a common practice 
here. [13 Nevertheless, to optimize the produc- 
tion process and to improve the homogeneity 
and strength of the fabricated workpiece one 
has to understand the complex rheological be- 
havior of the suspension and its relation to the 
microscopic structure. This knowledge in turn 
can be applied to soil mechanics. Shear thin- 
ning as observed in our simulations and exper- 
iments is an important mechanism for the dy- 
namics of landslides making them more danger- 
ous. 

In this paper we present our simulation re- 
sults of sheared suspensions of AI2O3 particles. 
The overall behavior is strongly determined by 
the effective interaction potential between the 
particles in the suspension. The potentials can 
be related to experimental conditions within 
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Debye-Hiickel theory, and thus we can compare 
our simulation results to experimental data. In 
contrast to our approach of a direct compari- 
son to experimental data, in the literature sim- 
ulation results are often compared to analytical 
calculations. 

Many different simulation methods have been 
developed and applied to colloidal suspen- 
sions: Stokesian dynamics fSD)[l8l IT9I |2Cl. 
accelerated Stokesian dynamics fASD)[2ll. |22|. 
pair dra g s imulations^."^. Brownian Dynam- 
ics (B DllM l25j . Lattice Boltzmann (LB)j2l 
I27L l28L l29l| an d Stochastic Rotation Dynamics 
( SRD) |30L l3lL l32l| . Due to the complex nature 
of the problem, all simulation methods have to 
simplify in some point. Either Brownian mo- 
tion is neglected or hydrodynamic interactions 
are included on a very simplified level. In many 
cases simulations are done without a quantita- 
tive comparison to experiments. In the present 
paper we combine Molecular Dynamics (MD) 
to simulate the colloidal particles, SRD for the 
description of the fluid, and a charge regulation 
model which provides us with realistic parame- 
ters for the Derjaguin Landau Vervey Overbeek 
(DLVO) potentials jH Hi] in the MD simula- 
tion. We include long range hydrodynamic in- 
teractions on a coarse grained level in the SRD 
part and we only include DLVO pair-potentials 
in the MD part. No electrostatic many body 
interactions or electrodynamic interactions are 
considered and modifications of the pair poten- 
tials due to locally increased colloid concentra- 
tions are neglected, too. However, many nu- 
merical investigations are based on much sim- 
pler models than ours. To our opinion our 
model covers the main properties quite well. 

Our paper is organized as follows: first we 
shortly describe our MD implementation fol- 
lowed by a short sketch of the SRD simula- 
tion method, and a description of how we have 
implemented our shear cell. The simulation 
method is described in detail in Ref . 3(J ■ Then 
we describe the so called 2pK charge regulation 
model which relates our simulation parameters 
with the pK- value and the ionic strength / ad- 
justed in the experiment. A short description 
of the simulation setup and of the experiments 
carried out follows. After that we present our 
simulation results and compare them to the ex- 
perimental data. Finally a summary is given. 



II. MOLECULAR DYNAMICS 

We study colloidal particles, composing the 
solid fraction, suspended in a fluid solvent. The 
colloidal particles are simulated with Molecular 



Dynamics (MD), whereas the solvent is mod- 
eled with Stochastic Rotation Dynamics (SRD) 
as described in Sec. lIIII 

In the MD part of our simulation we in- 
clude effective electrostatic interactions and van 
der Waals attraction, a lubrication force and 
Hertzian contact forces. The electrostatic and 
van der Waals potential are usually referred to 

as DLVO potentialsHHmHIiHlii, whidl 
capture the static properties of colloidal parti- 
cles in aqueous suspensions. The first compo- 
nent is the screened Coulomb term 
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where d denotes the particle diameter and r is 
the distance between the particle centers, e is 
the elementary charge, T the temperature, fee 
the Boltzmann constant, and z is the valency 
of the ions of added salt. Within DLVO the- 
ory one assumes linear screening, mainly by one 
species of ions with valency ±z (e.g. z = +1 for 
NHj"). The first fraction in Eq.^is a correction 
to the original DLVO potential, which takes the 
surface curvature into account and is valid for 
spherical particles|37|. 

The effective surface potential £ is the elec- 
trostatic potential at the border between the 
diffuse layer and the compact layer, it may 
therefore be identified with the ^-potential. It 
includes the effect of the bare charge of the col- 
loidal particle itself, as well as the charge of the 
ions in the Stern layer, where the ions are bound 
permanently to the colloidal particle. In other 
words, DLVO theory uses a renormalized sur- 
face charge, which we determine by the model 
described in Sec. llVI 

k is the inverse Debye length defined by 
k 2 = 8tt£bI, with the ionic strength /. The 

a 2 

Bjerrum length £g := 4 ^ £ e measures the dis- 
tance at which the electrostatic interaction of 
two elementary charges amounts /3 _1 = k^T. 
Sq is the permittivity of the vacuum, e r the rel- 
ative dielectric constant of the solvent (we use 
81 for water, i.e., Ib — 7 A for room tempera- 
ture). 

The Coulomb term of the DLVO potential 
competes with the attractive van der Waals 
term 
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A H = 4.76 • 1CT 20 J is the Hamaker constant plf 
which involves the polarizability of the parti- 
cles. The singularity of Vyaw for touching par- 
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tides is removed and the primary minimum is 
modeled by a parabola as described in Ref. [s3| ■ 
Long range hydrodynamic interactions are 
taken into account in the simulation for the 
fluid as described below. This can only repro- 
duce interactions correctly down to a certain 
length scale. On shorter distances, a lubrica- 
tion force has to be introduced explicitly in the 
MD simulation. The most dominant mode, the 
so-called squeezing mode, is an additional force 

Fl „ b =_ (w)f ^(!) 2 (3, 

between two particles with radius R and rel- 
ative velocity v ro i. 77 is the dynamic viscosity 
of the fluid. In contrast to the DLVO poten- 
tials the lubrication force is a dissipative force. 
When two particles approach each other very 
closely, this force becomes very large. To en- 
sure numerical stability of the simulation, one 
has to limit Fi u b- We chose a maximum force 
at a certain gap width r sc and shift the force so 
that the maximum force cannot be exceeded: 
Instead of calculating Fi u b(r) we take the value 
for Fi u b(r + r sc ). In addition, since the force 
decays for large particle distances, we can in- 
troduce a large cutoff radius r\ c for which we 
assume Fi u b(r) = if r — d > r\ c . As the in- 
tention of Fiub is to correct the finite resolution 
of the fluid simulation, r sc and r\ c have to be 
adjusted in a way that the dynamic properties, 
i.e., the viscosity of a simulated particle sus- 
pension with weak DLVO interactions fits the 
measurements. It turns out that r sc = and 
r \c = §rf work best. Our approach for Fi u b is 
similar to the one often used in lattice Boltz- 
mann simulations [2^. In contrast to LaddpCjj 
we have chosen to use two cut-off radii to be 
able to treat small and large gaps separately. 
There are different approaches, e.g., for Stoke- 
sian dynamics 01, where the force field is ex- 
panded to a multi pole series and the far field 
part is subtracted afterwards. 

Finally we use a Hertz force described by the 
potential 

Vfecrtz = K(d - rf' 2 if r<d, (4) 

where K is the constant which describes the 
elasticity of the particles in the simulation. The 
Hertz force avoids that the particles penetrate 
each other. It also contains a damping term in 
normal direction, 

F D am P = -(v re i,r)f/3Wd- r, (5) 
with a damping constant /3d. 



Since in this work no stress perpendicular to 
the shear direction is applied, the tangential 
forces at the particle surface are not of essential 
importance. To verify this, we have increased 
the spacial resolution of the fluid simulation, in- 
cluded tangential forces on the particles and al- 
lowed particle rotations. Even though the com- 
putational effort was considerably larger and 
one could expect that more effects on the length 
scale below the particle diameter could be cov- 
ered. However one could observe only a change 
of some percent in the viscosity and in the ve- 
locity profile. Due to the DLVO potential and 
the lubrication force the particles very rarely 
get into contact as long as no confining stress is 
applied. The only case, in which particles really 
touch each other, would be if the (^-potential 
is close to zero at a certain pK value. This 
pK value is called "isoelectric point" . It de- 
pends on the material of the suspended parti- 
cles and on the solvent. For our system it is at 
pH = 8.7|33. I n experiments close to the iso- 
electric point a solid fraction immediately floc- 
culates out and sediments. In the simulation 
one ends up with only one big cluster in the 
simulation volume, which corresponds to a part 
of a floe seen in the experiment. 

For this study we do not apply tangential 
forces and thus, having only central forces, we 
could neglect rotation of the particles. This re- 
duces the computational effort substantially. 



III. STOCHASTIC ROTATION 
DYNAMICS (SRD): SIMULATION OF 
THE FLUID 

The stochastic rotation dynamics method 
(SRD) was first introduced by Malevanets and 
Kapral [38l l39| . The method is also known 
as "real-coded lattice gas" |31| or as "multi- 
particle-collision dynamics" (MPCD) and has 
been successfully applied to simulate many im- 
portant systems like complex fluids containing 
polymers pol l4l|. vesicles in flow|42j. dynam- 
ics of chemical reactions 43j . The method is a 
promising tool for a coarse-grained description 
of a fluctuating solvent, e.g. in Ref. results 
of simulations of a flow around a cylinder are 
presented, or in Ref. [3^ sedimentation of a par- 
ticle suspension is studied. 

The method is based on so-called fluid par- 
ticles with continuous positions and velocities. 
Each time step is composed of two simple steps: 
One streaming step and one interaction step. 
In the streaming step the positions of the fluid 
particles are updated as in the Euler integra- 
tion scheme known from Molecular Dynamics 
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simulations, 

r i (t + T)=r i (t)+TV i (t), (6) 

where (t) denotes the position of the particle 
i at time t, Vi(t) its velocity at time t and r 
is the time step used for the SRD simulation. 
After updating the positions of all fluid parti- 
cles they interact collectively in an interaction 
step which is constructed to preserve momen- 
tum, energy and particle number. The fluid 
particles are sorted into cubic cells of a regular 
lattice and only the particles within the same 
cell are involved in the interaction step. First, 

their mean velocity Uj(i') = N l / t ,^ ^ v i(*) 

is calculated, where Uj(i') denotes the mean ve- 
locity of cell j containing Nj(t') fluid particles 
at time t' = t + r. Then, the velocities of each 
fluid particle in cell j are updated as: 

v i (t + r) = u J (t')+n i (t , )-[v i («)-u 3 - (*')]• (7) 

flj(t') is a rotation matrix, which is indepen- 
dently chosen randomly for each time step and 
each cell. We use rotations about one of the co- 
ordinate axes by an angle ±a, with a fixed |45|. 
The coordinate axis as well as the sign of the 
rotation are chosen at random, resulting in 6 
possible rotation matrices. The mean velocity 
Uj (t) in the cell j can be seen as streaming ve- 
locity of the fluid at the position of the cell j at 
the time t, whereas the difference [vj(t) — Uj(i')] 
entering the interaction step can be interpreted 
as a contribution to the thermal fluctuations. 
Thus, to calculate the local temperature in the 
cell under consideration one has to sum over 
the squares of this expression. 

The method just described is able to repro- 
duce hydrodynamics and thermal fluctuations. 
To couple the colloidal particles to the stream- 
ing field of the solvent, we use "Coupling II" 
of Ref. JsOl : we modify the rotation step of the 
original SRD algorithm slightly. The colloidal 
particles are sorted into the SRD cells as well 
and their velocity enters into the calculation of 
the mean velocity Uj (t) in cell j. Since the mass 
of the fluid particles is much smaller (in our case 
it is 250 times smaller) than the mass of the col- 
loidal particles, we have to use the mass of each 
particle — colloidal or fluid particle — as a weight 
factor when calculating the mean velocity 

%(*0 = E Vi(*)mi,(8) 

N 3 (t') 

with Mj(if) = ^2 m *> ( 9 ) 

i=l 



where we sum over all colloidal and fluid parti- 
cles in the cell, so that Nj(t') is the total num- 
ber of both particles together, is the mass 
of the particle with index i and therefore Mj(t') 
gives the total mass contained in cell j at the 
time t' = t + t. The update rule for the parti- 
cle velocities v,(t) and positions r^t+r), which 
we apply, is summarized in Eqns. © 0- This 
method to couple some embedded material to 
the SRD simulation is described for different 
applications in the literature 0,E3- 
This coupling method does not enforce no-slip 
boundaries on the particle surface, as in the 
method suggested by Inoue et al.|3l|. More- 
over, as very recently discussed by Padding and 
Louis purely radial interactions effectively 
introduce slip boundary conditions. Consider- 
ing the drag coefficient, a pre- factor changes 
and this could be corrected by assuming a dif- 
ferent hydrodynamic radius. We have checked 
the influence of hydrodynamic interactions by 
removing the fluid completely and by varying 
the resolution of the SRD simulation. Also 
two different coupling methods as described in 
Ref. have been applied. Without fluid, the 
achieved shear rate as well as the viscosity dif- 
fered strongly, whereas the difference between 
the two coupling methods was in the order of 
some per cent only. Therefore, we need hydro- 
dynamics to some extend, but we have chosen 
the coupling method with less computational 
effort. Very recently, Yamamoto et al. have 
shown that for colloidal gelation hydrodynamic 
interactions are of minor importance 49] for 3D 
systems, but in contrast to our work, they fo- 
cus on the static properties of a colloidal system 
quenched to zero temperature. 

In Ref. [3(J we have described a simple 
method to introduce shear at the fluid bound- 
ary by adding a velocity offset to all fluid par- 
ticles reflected at the shear plane. From a con- 
stant velocity offset Av one can calculate the 
mean shear force 




where L denotes the average number of fluid 
particles crossing through the shear plane in 
one time step and (. . .) stands for a time av- 
erage. L can be expressed by the mean free 
path and the number density of fluid particles. 
This would be a force driven shear, where one 
has only indirect control on the shear rate 7 or 
the shear velocity vg respectively. Therefore, 
we modify the mean velocity Uj(i') in the cells 
close to the shear plane by changing the veloc- 
ity of each fluid particle as well as the veloc- 
ity of the colloidal particles contained in that 



5 



specific cell by the difference vs — Uj(t'). By 
construction the mean velocity in these cells is 
equal to the shear velocity vg after that step. 
At the wall itself we implement full slip bound- 
ary conditions for the fluid and for the colloidal 
particles. The boundary in the direction of the 
shear profile (direction of the velocity gradient) 
is chosen to be non-periodic. By doing so, we 
can also observe phenomena like wall-slip, non- 
linear velocity profiles or density profiles in our 
shear cell (see Sec. I VII C|l . In the case of a non- 
linear velocity profile the viscosity is not well 
defined. We extract the central region of the 
profile where it is in first approximation lin- 
ear and estimate there an averaged viscosity. 
This is the ratio of the velocity gradient and 
the shear force which can be calculated in anal- 
ogy to Eq. I|l 0(1 by carrying out the sum over all 
velocity changes made. The region where we es- 
timate the velocity gradient is half the system 
size. 

We have tested a number of boundary con- 
ditions and different ways to impose shear, but 
the method just described turned out to work 
best. No-slip boundaries at a top and bottom 
plane seemed to work for high volume frac- 
tions and unless the potentials get attractive. 
As soon as (only slight!) cluster formation 
sets in the particles concentrate in the center 
of the system and lose contact to the sheared 
walls. Shearing only the fluid and not the col- 
loidal particles works always, but the result- 
ing viscosity is much too small. In fact, what 
one measures is the flow of the fluid stream- 
ing around the particles like a flow through a 
porous medium. The next point is how to deter- 
mine the shear force and the velocity gradient 
we need for the calculation of the shear viscos- 
ity. The force is always related to any velocity 
changes made in the system and its calculation 
is straightforward in most cases. 
The imposed velocity difference divided by the 
system size perpendicular to the shear plane 
would give an averaged gradient. For clustered 
systems not even the shape of the shear viscos- 
ity against shear rate was comparable to the 
measurements. If the velocity gradient changes 
within the system (compare Fig.^J, we have 
to take care that we measure the viscosity in 
the bulk, i.e., that we take the velocity gra- 
dient there. At least if the particles are not 
too strongly clustered, the slope of the plateau 
in the center of the system can be taken as a 
"good" velocity gradient. We use this velocity 
gradient as achieved shear rate as mentioned 
above. With this scheme for strongly attrac- 
tive forces the obtained viscosity 77(7) for the 
simulation stays in the vicinity of the measured 



curve, whereas for other methods we tried out 
the points of the simulation usually ended up 
far off the measured curve. 
Fully periodic boundary conditions for sheared 
systems, known as Lees Edwards conditions 
would be a good choice for stable suspensions. 
As soon as clusters are formed, the velocity 
profile becomes non linear, as discussed above. 
But, additionally the location of the cluster, 
i.e., the position of the plateau, is not fixed any- 
more to the center of the system, which makes 
it more difficult to extract the correct veloc- 
ity gradient. In addition, the shear force would 
be determined from the velocity changes of the 
particles passing around the periodic bound- 
aries. If the cluster by chance stays in the cen- 
ter of the system, again only the fluid would 
be sheared and only indirectly, transmitted by 
the fluid, the force would be exerted on the 
particles, as if with closed boundaries only the 
walls would move and no sheared regions close 
to the wall were implemented. Together with 
the periodic boundaries this would lead to large 
fluctuations of the shear force, caused by the 
present position of the cluster. Furthermore, 
the boundary conditions would have to be con- 
sistent for the MD and for the SRD simula- 
tion. For the MD part it is important that the 
position, where a particle re-enters the system 
after passing around the periodic boundary, is 
shifted by 2tv$ with t being the continuously 
increasing simulation time and the shear ve- 
locity. Additionally this shift has to be wrapped 
around the periodic boundaries in shear direc- 
tion. If we do the same with the fluid particles 
the shift could be any value, not necessarily an 
integer multiple of the fluid box size. What we 
want to point out, without any further restric- 
tions, the grid in the SRD rotation step would 
not anymore be regular in this plane, which in 
addition is the plane where one measures the 
shear force. To overcome this problem, one can 
restrict the shear rate to values determined by 
the SRD grid size and the SRD time step, but 
the other difficulties mentioned before remain. 



IV. THE CHARGE REGULATION 
MODEL 

To determine the effective surface potential 
which enters the DLVO potential, we use the 
model described in the following. In reality, the 
surface charge is achieved by adsorption and 
desorption of charge determining ions leading 
to an electrostatic potential difference between 
surface and bulk which in turn influences ion 
adsorption. A full description of this regula- 
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tion of surface charges requires two parts: the 
first part describes the relation between surface 
charge density and surface potential due to the 
electrolytic environment, whereas the second 
part quantifies the ion adsorption depending on 
the surface concentration of charge determining 
ions. 

Concerning the first part, a relation between 
the surface charge density a and the surface po- 
tential £ of a charged spherical colloidal particle 
of radius R immersed in an electrolytic environ- 
ment of relative dielectric constant e r and ionic 
strength I is given within Debye-Hiickel theory 

mnii by 



c 



Rcr 



eos r (l + kR) 



(11) 



As mentioned above, we consider the Stern 
layer as a part of the surface charge and thus, 
we can identify the effective surface potential 
in DLVO theory with the ^-potential and we 
can thus skip a discussion of bare charge versus 
effective charge [5^ EH 0^1 ■ 

In the second part of our model, the adsorp- 
tion of charge determining ions on the surface of 
the colloidal particle is described by assuming 
that the only mechanism of adsorption is that 
of protons (H + ) on surface sites (— S). It turned 
out that this assumption leads to reasonable re- 
sults for surfaces made of AI2O3. Adsorption is 
described by the two chemical reactions 55] 



■S - 
-SH 



H+ 
H+ 



-SH, 



-SH+ 



with the two reaction constants 



K x := 
K 2 := 



[-S-][H+]exp(-/feC) 
[-SH] 

[-SH][H+]exp(-/3eC) 

[-SH+] 



(12) 
(13) 



(14) 
(15) 



In terms of the surface site concentrations, 
the total number of surface sites per area 
and the surface charge density are given by 
N s = [-S-] + [— SH] + [-SH+] and a = 
— e[— S~] + e[— SHj"], respectively. Defining 
pK x := -\og w {Ki) and pK 2 := -log 10 (Jf 2 ) 
yields the point of zero charge pH z , i.e., the 
pH value of vanishing surface charge, as pH z = 
h{pK\ +PK2). The surface site density Ns and 
the difference ApK := pK\ — pK 2 are treated 
as adjustable parameters. 

The above equations lead to the relation 



a 8 sinh(-07v — /3e£) 

eN s ~ l + (5cosh(7/-Ar-/?eC) 



(16) 



with the Nernst potential ipjy '■= ln(10)(pH z — 
pH) and 6:= 2- 10"^. 

Equations fTJ and H16[l can be solved self- 
consistently for £ as a function of pH. For 
our system of AI2O3 particles we find ApK = 
4.2 and N s = 0.22/nm 2 . With these values 
the measured ^-potential of 52 mV at pH = 
6, I = 0.01mol/l and up to 120 mV at 
pH = 4, I = 0.01mol/l can be reproduced 
best. For the experimental determination of 
the ^-potential electrophoretic (Delsa 440SX, 
Beckman-Coulter GmbH, Germany) and elec- 
trokinetic measurements (AcustoSixer lis, Col- 
loidal Dynamics Ind., USA) were performed. 
To calculate the £-potential Henry's theory^] 
was used. For details see Ref. m We have 
to admit that the relation between the directly 
measured quantities, e.g., electrophoretic mo- 
bility, and the ^-potential is subject of current 
research[53IH|5!|I3. 



V. SIMULATION SETUP 

In our simulation we try to model the ex- 
perimental system as good as possible. We 
start with spherical particles of diameter d = 
0.37 pm, the mean diameter of the particles 
used in the experiment. The simulation box 
is 48d = 17.76 pm long in x-direction, 24d = 
8.88 /xm in z-direction, and 12d = 4.44 /xm in 
y-direction. To achieve a volume fraction of 
usually $ = 35%, as in the experiment, we 
need to simulate simulate 9241 spheres. Our 
shear direction is the x-direction and the ve- 
locity gradient of the shear flow points in z- 
direction; in other words we shear the upper 
and lower xy-plane with respect to each other 
in x-direction. We use periodic boundaries in 
x- and y-direction and closed boundaries in 
z-direction for both, fluid and MD particles. 
The energy supplied by the shear force is dis- 
sipated by means of a Monte-Carlo-thermostat 
described in Ref. [3(J H3l- H acts on the fluid 
particles as well as on the MD particles and 
conserves the momentum in each SRD cell. 



VI. EXPERIMENTAL SETUP 

Experiments are carried out with high pu- 
rity (99,97%) a-Al 2 3 powder (RCHP DBM, 
Baikowski Malakoff Industries, Inc., USA). The 
mean particle diameter is 0.367/im (Coulter LS 
Particle Size Analyzer) and the size distribution 
is narrow (dio = 0.176/im, dgo — 0.664pm). 
The powder is suspended in bidistilled water 
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(Merck, Germany). The suspension is then dis- 
persed with alumina balls in a ceramic con- 
tainer for 24 h at a small rotational speed to 
keep the abrasion low. Subsequently, the sus- 
pension is degassed at 50 mbar under agitation. 
Then, in order to reduce the ionic strength 
to the desired degree, the suspension is pu- 
rified by the dialysis technique. In this way 
the majority of ions are removed and a back- 
ground electrolyte of a very low salt concen- 
tration (5 • 10 _4 mol/l) is obtained for suspen- 
sions of high solids loading. Starting from this 
master suspension, suspensions with increased 
ionic strength are obtained by adding differ- 
ent amounts of dry ammonium chloride NH4CI 
(Merck, Germany). The pH of the suspensions 
is adjusted to pH = 6 with 0.1 and lmol/1 hy- 
drochloric acid HC1 (Merck, Germany), if nec- 
essary. Thereby the ionic strength and pH are 
revised by use of a laboratory pH- and conduc- 
tivity meter (inoLab pH/Cond Level 2, WTW 
GmbH, Germany). The electrophoretic mo- 
bility of dilute suspensions is measured with 
a Coulter Delsa 440 SX. Irreversible aggrega- 
tion due to inhomogeneous salt concentration 
is not of importance here. If the ionic strength 
is strongly increased and after that, a second 
dialysis step is performed to remove the ions 
again, the original viscosity is restored. 

The ion concentrations of selected ions are 
measured before and after dialysis using Induc- 
tively Coupled Plasma-Optical Emission Spec- 
troscopy (ICP-OES, Model JY 70 plus, France). 
The suspensions are characterized using a Vis- 
colab LC10 rheometer (Viscolab LM rheome- 
ter with control unit Viscolab LC10, Physica, 
Germany) with a cup and bob or a double gap 
geometry. The measurements are either per- 
formed immediately after suspension prepara- 
tion, or they are stored on a roller bank to avoid 
sedimentation. Sedimentation during the ex- 
periment can be excluded, since it takes much 
longer than the whole experiment, and the 
shear forces are much larger than gravity. Shear 
rate controlled experiments are performed at 
a constant temperature of 20° C. The suspen- 
sions are sheared at a constant shear rate of 
7 = 300/s (Pe = 8.8) before starting the ac- 
tual ramp measurement. In the experiments 
the shear rate is increased up to 7 = 4000 / s 
(Pe — 117) and decreased again to zero. In 
this paper, whenever referring to a shear rate 
we also specify the Peclet number 

Pe = 6irriR 3 j/k B T, (17) 

to make it easier for the reader to compare our 
results to other data. When the suspensions are 
pre sheared an occurring discrepancy between 



100 







+^ clustered 




V B 








- suspended \ 


/suspended 

/ s' 


repulsive x \ 

II I N I I I I 


/ . 'repulsive 
li I 



3 4 5 6 7 8 9 10 11 12 13 14 
pH value 

FIG. 1: Schematic stability diagram for volume 
fraction $ = 35% in terms of pH-value and ionic 
strength involving three different microstructures: 
A clustering regime due to van der Waals attrac- 
tion, stable suspensions where the charge of the col- 
loidal particles prevents clustering, and a repulsive 
structure for further increased electrostatic repul- 
sion. This work concentrates on state A (pH — 6, 
/ = 3mmol/l) in the suspended region, state B 
(pH = 6, I = 7mmol/l) close to the border but al- 
ready in the clustered region, and state C (pH = 6, 
/ = 25mmol/l) in the clustered region. The bor- 
ders are not sharp transitions, but notable in a 
change of the shear viscosity. 



the measured viscosity in the increasing ramp 
and the decreasing one can be minimized. A 
detailed description of the experiments will be 
published elsewhere [Til lolf . 



VII. RESULTS 

A. Stability Diagram 

Depending on the experimental conditions, 
one can obtain three different microstructures: 
A clustered region, a suspended region, and 
a repulsive structure. The charge regulation 
model allows us to quantitatively relate the 
interaction potentials to certain experimental 
conditions. A schematic picture of the stabil- 
ity diagram is shown in Fig.^ Close to the 
isoelectric point (pH = 8.7), the particles form 
clusters for all ionic strengths since they are 
not charged. At lower or higher pH values one 
can prepare a stable suspension for low ionic 
strengths because of the charge, which is carried 
by the colloidal particles. At even more extreme 
pH values, one can obtain a repulsive structure 
due to very strong electrostatic potentials (up 
to C = 170 mV for pH = 4 and 1=1 mmol/1, 
according to our model). The repulsive struc- 
ture is characterized by an increased shear vis- 
cosity. In the following we focus on three states: 
State A (pH = 6, I — 3 mmol/1) is in the sus- 
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FIG. 2: (color online) Images of four different cases. 
For better visibility we have chosen smaller systems 
than we usually use for the calculation of the vis- 
cosity. The colors denote velocities: Dark particles 
are slow, bright ones move fast. The potentials do 
not correspond exactly to the cases A—C in Fig.0 
but they show qualitatively the differences between 
the different states: 

a) Suspension like in state A, at low shear rates. 

b) Layer formation, which occurs in the suspension 
(state A) at high shear rates and in the repulsive 
regime already at moderate shear rates. 

c) Strong clustering, like in state C, so that the sin- 
gle cluster in the simulation is formed. 

d) Weak clustering close to the border like in state 
B, where the cluster can be broken into pieces, 
which follow the flow of the fluid (plug flow). 



pended region, state B (pH = 6, I = 7mmol/l) 
is a point already in the clustered region but 
still close to the border, and state C (pH = 6, 
/ = 25mmol/l) is located well in the clustered 
region. 

Some typical examples for the different mi- 
crostructures are shown in Fig.^,)— d). These 
examples are meant to be illustrative only and 
do not correspond exactly to the cases A-C 
in Fig.^ denoted by uppercase letters. In the 
suspended case (a), the particles are mainly 
coupled by hydrodynamic interactions. One 
can find a linear velocity profile and a slight 
shear thinning. If one increases the shear rate 
7 > 500/s (Pe > 15) , the particles arrange in 
layers. The same can be observed if the Debye- 
screening length of the electrostatic potential 
is increased (b), which means that the solvent 
contains less ions (/ < 0.3 mmol/1) to screen the 
particle charges. On the other hand, if one in- 
creases the salt concentration, electrostatic re- 
pulsion is screened even more and attractive 
van der Waals interaction becomes dominant 
(/ > 4 mmol/1). Then the particles form clus- 
ters (c), and viscosity rises. A special case, 



called "plug flow", can be observed for high 
shear rates, where it is possible to tear the clus- 
ters apart and smaller parts of them follow with 
the flow of the solvent (d). This happens in our 
simulations for / = 25 mmol/1 (state C) at a 
shear rate of 7 > 500/s (Pe > 15). However, as 
long as there are only one or two big clusters in 
the system, it is too small to expect quantita- 
tive agreement with experiments. In these cases 
we have to focus on state B (I = 7 mmol/1) 
close to the border. 

In our simulations we restrict ourselves to the 
region around pH = 6 where we find the border 
between the suspended region and the clustered 
regime at about I = 4 mmol/1 in the simula- 
tions as well as in the experiments. Also the 
shear rate dependence of the viscosity is com- 
parable in simulations and experiments as dis- 
cussed in Sec. lVII CI 

In Ref. [62( a qualitative stability diagram 
similar to Fig.^ has been shown. The bor- 
ders there are shifted, since they depend on the 
threshold value for which one defines that the 
viscosity has increased. Correspondingly, if one 
is less sensitive on the viscosity increase, one 
would still consider the system to be suspended, 
if only weak cluster formation takes place. 



B. Total Energy 

In our simulations we calculate the total en- 
ergy, because it can be used as a tool to check 
if the response of the simulation to the vari- 
ation of any parameter is consistent with the 
expectations, e.g., a decrease of the surface 
charge on the colloidal particles should cause 
the secondary minimum of the DLVO potential 
to become deeper and thus decrease the total 
energy-but if the total energy increases, this 
can be an indication for numerical instabilities. 

The total energy comprises the kinetic energy 
of both fluid and colloidal particles, including 
thermal motion on the microscopic level, as well 
as the potential energy due to Coulomb repul- 
sion, van der Waals attraction, and Hertz con- 
tact forces. Our simulations are carried out at 
room temperature (T = 295 K) and constant 
volume fraction. Supposing a linear velocity 
profile, the kinetic energy increases quadrati- 
cally with the shear rate 7. This can be ob- 
served if the electrostatic repulsion is, on the 
one hand, strong enough to prevent cluster for- 
mation due to van der Waals attraction and, on 
the other hand, weak enough, so that the col- 
loidal particles can move relatively freely with- 
out undergoing a glass transition or crystalliza- 
tion. 
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FIG. 3: Total energy depending on the shear rate 
7 for the states A (I — 3mmol/l) and C (I = 
25 mmol/1) of Fig.Q In state A the system is a sta- 
ble suspension, in state C cluster formation reduces 
the total energy at low shear rates. At 7 = 500/s 
(Pe = 15), the cluster can be broken up into two 
parts moving in opposite directions. The two solid 
bodies have a larger kinetic energy than the suspen- 
sion with a linear velocity profile. This explains the 
crossover of the two curves. For even higher 7, the 
clusters are broken up in more pieces leading ulti- 
mately to the same structure as for the suspended 
state A. The energy axis has been scaled by the 
total number of particles (fluid particles plus col- 
loidal particles) and plotted in units of ksT. The 
solid line is the analytical solution (Eg . 11811 for a 
linear velocity profile, the dashed lines are a guide 
to the eye. 



If the interactions are strongly repulsive, i.e., 
in the case of very low salt concentration, where 
the Debye-screening length is large, one can see 
an extra contribution of the electrostatic repul- 
sion to the total energy. If the volume frac- 
tion is low, the particles can still find a con- 
figuration, in which the mean nearest neighbor 
distance is larger than the interaction range of 
the repulsion. But, if the volume fraction is in- 
creased, the particles have to be packed closer, 
which leads to a constant positive offset to the 
total energy. It only depends on the poten- 
tials and on the volume fraction, but not on 
the shear rate. 

In a similar way as for repulsive interactions, 
one can understand the negative energy con- 
tribution in the case of high salt concentra- 
tions: The DLVO potentials contain a min- 
imum where attractive Van der Waals inter- 
action is stronger than electrostatic repulsion. 
Then the particles form clusters and "try" to 
minimize their energy. In Fig. [3] for small shear 
rates the values for the energy in the clustered 
case of state C (I = 25 mmol/1) is lower than for 
the suspended case of state A (I = 3 mmol/1). 
We have plotted the total energy divided by the 



total number of all particles (fluid particles plus 
colloidal particles) in units of /cbT. For 7 — > 
the energy per particle approaches ^k#T as one 
would expect in 3D. The solid line in Fig. [3] is 
the sum of §/cbT per particle with the kinetic 
energy of a fluid with a linear velocity profile: 

E tot = ^k B TN tot + ^Vgj 2 Ll (18) 

where N to t is the total number of both, fluid 
and colloidal particles, V denotes the volume 
of the simulated system, g the averaged mass 
density of the suspension, and L z is the exten- 
sion in z-direction (perpendicular to the shear 
plane). State A coincides very well with this 
curve. 

For state C the behavior is shear rate depen- 
dent: In contrast to the repulsive case, clusters 
can be broken up. This happens at a shear 
rate 7 = 500/s (Pe = 15) (see Fig.OJl where 
one obtains two clusters moving in opposite di- 
rections. Since in this case the resistance of 
the system decreases, the velocity of the two 
clusters becomes larger. Since both clusters are 
moved as a whole, their energy even becomes 
larger than in the suspended case. If one fur- 
ther increases the shear rate, no (big) clusters 
can form anymore and the energies for both salt 
concentrations are nearly the same, and corre- 
spond to the kinetic energy of a suspension with 
a nearly linear velocity profile. For 7 = 1000/s 
(Pe — 29) state C coincides with the analytic 
curve. 



C. Shear Profile and Shear Viscosity 

In each of the three regimes a typical ve- 
locity profile of the shear flow occurs. For 
the suspended microstructure one finds a lin- 
ear velocity profile (Fig.Qti)) with nearly New- 
tonian flow. The particles are distributed ho- 
mogeneously, thus the density profile is struc- 
tureless (Fig.EJi)). The motion of the par- 
ticles is only weakly coupled by the hydro- 
dynamic forces. At high enough shear rates 
(7 > 500/s (Pe > 15)) the particles arrange 
in layers parallel to the shear plane, as can be 
seen in the density profile Fig.^b), too. This 
arrangement minimizes collisions between the 
particles. As a result, the shear viscosity de- 
scents as shown in Fig-EI which we discuss more 
in detail below. Shear induced layer forma- 
tion has been re port ed in literature for different 
experiments |63t l64l l6 5| . |6 6| . |67l and Stokesian 
dynamics simulations [1^, log . For low shear 
rates Brownian motion disturbs the layers or 
prevents their formation. As shown in Ref. [l9| , 
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FIG. 4: Profiles of tangential velocity component 
(v x ) in normal direction (2): 

a) Linear profile in the suspended regime, state A 
of Fig.[U(/ = 3mmol/l) at 7 = 500/s (Pe = 15)). 

b) Cluster formation in state C (I = 25mmol/l) at 
7 = 100/s (Pe — 2.9). In principle one could de- 
termine the viscosity of one single cluster from the 
central plateau, but this is not the viscosity found 
in experiments. There, one measures the viscosity 
of a paste consisting of many of these clusters. 

c) Same as case b) but with higher shear rate 
(7 = 500/s Pe = 15). Hydrodynamic forces are 
large enough to break the cluster into two pieces. 
The velocity axis is scaled with the shear velocity 
vs for better comparability. 



hydrodynamic forces can destroy them as well, 
if the shear rate is high enough. In our simula- 
tions we do not reach these conditions. In the 
simulations shear rates up to 2000/s (Pe = 59) 
can be realized before limitations of the simu- 
lation method influence the results. The incre- 
ment of the shear angle in one SRD time step 
7 = ji~sRD amounts about ir/4 then, i.e., the 
offset in x-direction between two neighboring 
layers of SRD cells in z-direction amounts one 
cell per SRD time step. 

Furthermore, the volume fraction and the in- 
teraction range of the electrostatic repulsion, 
or the ionic strength respectively, influence the 
layer formation: In the repulsive regime the lay- 
ers are more clearly already at moderate shear 
rates. It can be excluded that the effect is 
purely a finite size effect, since for un-sheared 
suspensions no layers can be observed, at least 
some particle diameters from the walls. In the 
repulsive regime the particles try to optimize 
their local structure, but a long range order as 
in the case of a sheared system can not be seen. 

In the clustered regime, the clusters move in 
the fluid as a whole. They are deformed, but 
since the inter-particle forces are stronger than 
the hydrodynamic forces, the cluster moves 
more like a solid body than like a fluid. Often 
there is one big cluster that spans the whole 



system. The density profile (FigJ^t)) increases 
in the central region and decays at the regions 
close to the border, since particles from there 
join the central cluster. When averaging the 
velocity profile in the shear flow, one finds a 
very small velocity gradient in the center of the 
shear cell and fast moving particles close to the 
wall, where the shear is imposed (Fig.^J))). The 
velocity profile is non-linear on the length scale 
of the simulations. In the experiment the phys- 
ical dimensions are much larger and therefore 
the velocity profile can become approximately 
linear again if the system consists of many large 
clusters. However, due to the computational 
effort in simulations it is today impossible to 
measure the shear viscosity for these strongly 
inhomogeneous systems. We have scaled our 
system by a factor of 2 in x and z-direction 
(keeping the volume fraction <3> = 35% con- 
stant), but we still observe one big cluster after 
some hundreds of SRD time steps, i.e., finite 
size effects are still present in our simulations. 

Closer to the border clusters can then be bro- 
ken up into small pieces by the hydrodynamic 
forces at least for high shear rates. In state 
C of Fig.H this happens for the first time at 
7 = 500/s (Pe = 15), so that one can find 
two clusters in the system moving in opposite 
directions. The velocity profile of this case is 
shown in Fig.QJ;). For even higher shear rates 
or closer to the border (e.g. state J5), the clus- 
ters are broken into smaller pieces. Then, they 
move in the shear flow with an approximately 
linear velocity profile. Due to van der Waals at- 
traction the system resists with stronger shear 
forces and the viscosity is higher than in the 
suspended case (Fig.^. 

In Fig. [S] the simulation results are shown to- 
gether with the experimental results, both for 
the two cases of a slightly clustered system in 
state B (I — 7mmol/l) and a suspension (state 
A, I = 3mmol/l). For the suspension (state 
^4) the viscosity decreases with the shear rate 
( "shear thinning" ) . The experimental data and 
the simulation are consistent within the accu- 
racy of our model. There are several reasons 
for which our model does not fit exactly the 
measurements: The most insecure factor which 
enters into the comparison is the measurement 
of the ^-potential. Starting from this point we 
set up our charge regulation model to extrap- 
olate to different salt concentrations, assum- 
ing two reactions being the only processes that 
determine the surface charge of the colloidal 
particles. Furthermore, we have monodisperse 
spherical particles, which is another simplifi- 
cation in our model. Then, the lubrication 
force as a correction for the finite resolution of 
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FIG. 5: Density profiles: a) Suspended case: State A in Fig.0(/ = 3mmol/l), at low shear rates (7 = 50/s 
(Pe = 1.5)). The density distribution is homogeneous. 

b) Shear induced layer formation: This is state A as in graph a) of this figure, but for a high shear rate 
(7= 1000/s (Pe = 29)). 

c) Strong attractive forces in state C (I — 25mmol/l): For low shear rates (7 = 50/s (Pe = 1.5)) only one 
central cluster is formed, which is deformed slowly. 
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FIG. 6: Comparison between simulation and exper- 
iment: viscosity in dependence of the shear rate, for 
the states A (I = 3 mmol/l) and B (I — 7 mmol/l) 
of Fig.Q Note: shear thinning is more pronounced 
for the slightly attractive interactions in state B 
than for the suspended state A. Lines denote ex- 
perimental data [rj|. points are results from our 
simulations. 



the fluid method can only recover to a certain 
degree the hydrodynamics on smaller length 
scales than the cell size of the fluid simulation, 
e.g., we have not implemented other modes of 
lubrication than the "squeezing mode" (Eq.[3J). 

However, we have done several tests where 
we have simulated systems with size and charge 
polydispersity in the order of magnitude corre- 
sponding to the experimental conditions. We 
have tested different sorts of boundary condi- 
tions, different ways to implement shear, and 
different coupling methods between fluid and 
particles. In most of the tests the achieved 
shear viscosity in the simulation did not change 
notably. To our experience the way shear is 
imposed and the particle size have the largest 
influence on the result. 

Finally, one has to keep in mind that the vis- 
cosity of the suspension can be varied by more 
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FIG. 7: Discrepancy between simulated viscosity 
and measurement for states A and B of Fig.0 
for different system sizes at low shear rates. The 
plot shows squared relative differences against z- 
extension of the simulation volume. The lines are 
a guide to the eye. 



than one order of magnitude, e.g., by changing 
the ionic strength. In this context the devi- 
ations between simulation and experiment are 
small. 

For the slightly clustered case (state B) an 
increase of the shear viscosity, compared to the 
suspended case, can be observed in the experi- 
ment as well as in the simulations. Shear thin- 
ning becomes more pronounced, because clus- 
ters are broken up, as mentioned above. How- 
ever, the shear rate dependence is stronger in 
the simulations than in the experiment. This 
can be the first indication of finite size effects. 

We have studied the dependence of the sim- 
ulated shear viscosity in dependence of the sys- 
tem size. The effect is most important for low 
shear rates and thus we carried out several sim- 
ulations for state A at 7 = 20/s (Pe = 0.6) 
and for state B at 7 = 50/s (Pe = 1.4) We 
have chosen these values because clustering is 
already too strong in state B at 7 = 20/s 
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(Pe = 0.6) to reasonably determine a viscos- 
ity and the system size dependence becomes 
too small for 7 = 50/s (Pe — 1.4) in state 
A. In Fig. [7| we plot the squared relative de- 
viation between simulation and measurements 
against the system size. We do not know, if 
the simulation results would exactly converge 
to the measured values if the simulated system 
is large enough. However, the figure shows the 
trend that the deviation becomes smaller for 
larger system sizes, but to reach in state B the 
same accuracy as in state A one would have 
at least to double the system size in each di- 
mension. It then takes approximately twice as 
long for the system to relax to a steady state, 
resulting in a factor of 16 in the computational 
effort. This gives a guess, that each single point 
of Fig.H3 would need approximately 3000 CPU 
hours. For smaller shear rates or even deeper 
in the clustered regime of the stability diagram, 
e.g., in state C (I = 25mmol/l), the finite size 
effects become more pronounced — ending up in 
the extreme case of only one big cluster exist- 
ing in the system. For simulations with good 
accuracy the effort again increases at least by 
the same factor. 

Unfortunately, exactly this would be the 
most interesting case with respect to soil me- 
chanics and to understand landslides, which 
was our initial motivation. Anyhow, if we com- 
pare state A and state B, shear thinning be- 
comes stronger with increasing ionic strength. 
In the experiments the effect becomes much 
stronger for larger ionic strengths (up to I = 
65mmol/l), where the viscosity for low shear 
rates is increased by more than a factor of ten. 
However, the fact that there is shear thinning 
and that it depends on the ionic strength is an 
interesting result, if we recall the fact that the 
lubrication force in our simulation can be inter- 
preted as a velocity dependent damping force, 
which becomes stronger for higher relative ve- 
locities. Therefore one would intuitively expect 
shear thickening 

Finally, the limitations of the DLVO the- 
ory have to be taken into account. DLVO 
potentials are derived for dilute suspensions 
and hence large particle distances. This is 
not fulfilled in our case-and even less inside 
the clusters. There are theoretical attempts 
that address the shortcomings of DLVO the- 
ory: Explicit simulation of micro- ions |69j. den- 
sity functional theory [toL FH F^ . response 
theory F3I Fil FBI Ft^ . Poisson Boltzmann cell 
model s I77I FM Fil. and full Poisson Boltzmann 
theory [80l l8ll l82l l83j |. but they have other 
disadvantages-most of them require a large 
computational effort. For Poisson Boltzmann 



cell models one assumes homogeneously dis- 
tributed colloidal particles, so that each of them 
can be regarded as a representative single par- 
ticle in a Wigner-Seitz cell. Additionally, de- 
pending on the level Poisson Boltzmann theory 
is included, a mysterious phase separation could 
be identified as an artifact of linearization[79|. 
Full Poisson Boltzmann theory would require 
the calculation of the local potential, not only 
as done in our charge regulation model in the 
beginning of the simulation, but for the whole 
simulation box and in each time step. This 
would in principle provide a better description 
of the real system, but the computational effort 
would be much larger making simulations of 
several thousands of particles impossible. The 
same applies to the approach of including the 
micro-ions explicitly in the simulation. One 
could obtain three (many) body interactions 
from full Poisson Boltzmann theory and try to 
include them as a lookup table in the simula- 
tion. However, one would have to decrease the 
system size to keep the computational effort 
affordable. For our simulations we need rela- 
tively simple pair potentials to keep the com- 
putational costs within a limit. Nevertheless, 
the overall behavior can be reproduced by the 
simulation on a semiquantitative level. The rea- 
son for that might be the fact that in some of 
the above mentioned theoretical attempts (den- 
sity functional theory and cell models) DLVO- 
like potentials are obtained with a renormalized 
charge and screening length. In our charge reg- 
ulation model we do nothing else than adjusting 
a renormalized charge to the measurements of 
the ^-potential. This may be a general expla- 
nation why DLVO potentials can often be used 
although the assumptions for DLVO theory are 
not fulfilledjH |H M, H HI . 

We have carried out simulations in the repul- 
sive region of the stability diagram as well. We 
find layers parallel to the shear plane in anal- 
ogy to FigJSJ}). In contrast to the suspended 
regime, in the repulsive regime the layer struc- 
ture is present — at least locally, but orientation- 
ally disordered — even if no shear is applied. If 
shear flow is present, the shear plane marks one 
orientation which the layer structure adopts. In 
some cases for very low ionic strengths one can 
observe shear bands so that the velocity gradi- 
ent and thus the viscosity as well vary strongly 
in the system. Again, in the experiment, phys- 
ical dimensions are much larger and on that 
length scale the velocity profile might be as- 
sumed to be linear when enough shear bands 
are in the system. The shear force and hence 
the viscosity rise with respect to the suspended 
regime, due to electrostatic repulsion. One can 
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FIG. 8: Viscosity versus volume fraction for the 
repulsive region (pH = 6 and I — 0.3mmol/l). The 
shear rate was 7 = 100/s (Pe = 2.9). The points 
are simulation results, the line is a guide to the eye. 



consider the particles together with the interac- 
tion range as soft spheres with an effective ra- 
dius of the interaction range of the electrostatic 
repulsion. This effective radius can be in our 
case about 25% larger than the particle radius. 
Therefore, a transition to a repulsive structure 
occurs in our systems already between 35% and 
40% volume fraction. Because of the smooth 
shape of the exponentially screened Coulomb 
potential it is not a sharp glass transition as for 
hard spheres, but smooth and shear rate de- 
pendent as well. In Fig.|H| we have shown the 
dependence of the viscosity on the volume frac- 
tion for pH — 6 and I = 0.3mmol/l. Starting 
at $ = 0.3 the shear viscosity starts to increase 
and reaches a value one decade larger beyond 
$ = 0.4. 



VIII. SUMMARY AND OUTLOOK 

We have shown how to relate DLVO poten- 
tials to the conditions in a real aqueous suspen- 
sion of AI2O3 particles. The behavior of shear 
viscosity has been studied in experiments and in 
simulations. We have found shear thinning due 
to a layer formation on the microscopic scale in 
the case of a suspension. 

If a clustered system is sheared, clusters are 
broken up into pieces by the imposed shear, 
which leads to a stronger shear thinning than in 
the suspended case. Close to the border we are 
able to reproduce the measured shear viscosity 



in the simulation. 

Deep in the clustered regime we have found 
that our particles form one big cluster in the 
system which can be broken up by the hydro- 
dynamic forces of the shear flow. For strongly 
clustered systems at low shear rate, which 
would be the most interesting case for soil me- 
chanics, there are strong finite size effects. One 
attempt to address this problem is to increase 
the size of the simulated system. As we have 
shown with our work, the computational effort 
increases to an extend that a parallelized simu- 
lation code would be necessary. However, this 
might be not sufficient, since in both cases, for 
low shear rates and for strongly attractive inter- 
actions the finite size effects become stronger. 
Since in these simulations a considerable am- 
mount of the computing time is consumed by 
the particles inside a cluster, one could think of 
a more coarse-grained description of the clus- 
ters. Nevertheless, input data for such a model 
could be obtained using our present simulation 
code. Depending on the model one might need 
information about the shear resistance of a sin- 
gle cluster, depending on the cluster size and 
shape. This would be very difficult to measure, 
but it could be calculated in a small simulation 
using our combined MD and SRD algorithm. 
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